clear all
load results.mat

% wholemat = [k,c,inv,lambda,PR Penergy]
%             1 2 3      4    5    6
% renewmat = [kr cr ir lambdar PRr P_Er Hr jr]
%              1  2  3  4       5    6  7  8
% fossilmat = [kf cf invf lambdaf PRf P_Ef Sf Nf nf gf sigmaf];
%              1   2   3     4     5   6    7  8  9 10   11


%% smooth N and n
smN0 = csaps(1:Tmat0(2)+1,fossilmat0(:,8),0.1,50:Tmat0(2)-2);
smN1 = csaps(1:Tmat1(2)+1,fossilmat1(:,8),0.1,50:Tmat1(2)-2);
smN2 = csaps(1:Tmat2(2)+1,fossilmat2(:,8),0.1,50:Tmat2(2)-2);
smN3 = csaps(1:Tmat3(2)+1,fossilmat3(:,8),0.1,50:Tmat3(2)-2);
plot([smN0',fossilmat0(50:Tmat0(2)-2,8)])

csmthn0 = smN0(2:end)-smN0(1:end-1);
csmthn1 = smN1(2:end)-smN1(1:end-1);
csmthn2 = smN2(2:end)-smN2(1:end-1);
csmthn3 = smN3(2:end)-smN3(1:end-1);
plot([csmthn0',fossilmat0(51:Tmat0(2)-2,9)]);

%% smooth k and i
smk0 = csaps(1:Tmat0(2)+1,fossilmat0(:,1),0.5,50:Tmat0(2)-2);
smk1 = csaps(1:Tmat1(2)+1,fossilmat1(:,1),0.5,50:Tmat1(2)-2);
smk2 = csaps(1:Tmat2(2)+1,fossilmat2(:,1),0.5,50:Tmat2(2)-2);
smk3 = csaps(1:Tmat3(2)+1,fossilmat3(:,1),0.5,50:Tmat3(2)-2);

csmthi0 = smk0(2:end)-(1-0.04)*smk0(1:end-1);
csmthi1 = smk1(2:end)-(1-0.04)*smk1(1:end-1);
csmthi2 = smk2(2:end)-(1-0.04)*smk2(1:end-1);
csmthi3 = smk3(2:end)-(1-0.04)*smk3(1:end-1);
plot([csmthi0',fossilmat0(51:Tmat0(2)-2,3)]);

%% smooth g
smg0 = (fossilmat0(50:end-3,10)+fossilmat0(51:Tmat0(2)-1,10))/2;
smg1 = (fossilmat1(50:end-3,10)+fossilmat1(51:Tmat1(2)-1,10))/2;
smg2 = (fossilmat2(50:end-3,10)+fossilmat2(51:Tmat2(2)-1,10))/2;
smg3 = (fossilmat3(50:end-3,10)+fossilmat3(51:Tmat3(2)-1,10))/2;
plot([smg1,fossilmat1(50:end-3,10)])

%% smooth PR and P_Er

PR0 = csaps(1:Tmat0(2)+1,fossilmat0(:,5),0.1,50:Tmat0(2));
PR1 = csaps(1:Tmat1(2)+1,fossilmat1(:,5),0.1,50:Tmat1(2));
PR2 = csaps(1:Tmat2(2)+1,fossilmat2(:,5),0.1,50:Tmat2(2));
PR3 = csaps(1:Tmat3(2)+1,fossilmat3(:,5),0.1,50:Tmat3(2));
plot([PR0',fossilmat0(50:Tmat0(2),5)]);

PE0 = csaps(1:Tmat0(2)+1,fossilmat0(:,6),0.1,50:Tmat0(2));
PE1 = csaps(1:Tmat1(2)+1,fossilmat1(:,6),0.1,50:Tmat1(2));
PE2 = csaps(1:Tmat2(2)+1,fossilmat2(:,6),0.1,50:Tmat2(2));
PE3 = csaps(1:Tmat3(2)+1,fossilmat3(:,6),0.1,50:Tmat3(2));
plot([PE0',fossilmat0(50:Tmat0(2),6)]);


%% replace data
fossilmat0(50:Tmat0(2)-2,[1 8 10]) = [smk0' smN0' smg0];
fossilmat1(50:Tmat1(2)-2,[1 8 10]) = [smk1' smN1' smg1];
fossilmat2(50:Tmat2(2)-2,[1 8 10]) = [smk2' smN2' smg2];
fossilmat3(50:Tmat3(2)-2,[1 8 10]) = [smk3' smN3' smg3];

fossilmat0(51:Tmat0(2)-2,[3 9]) = [csmthi0' csmthn0'];
fossilmat1(51:Tmat1(2)-2,[3 9]) = [csmthi1' csmthn1'];
fossilmat2(51:Tmat2(2)-2,[3 9]) = [csmthi2' csmthn2'];
fossilmat3(51:Tmat3(2)-2,[3 9]) = [csmthi3' csmthn3'];

fossilmat0(50:Tmat0(2),[5 6]) = [PR0' PE0'];
fossilmat1(50:Tmat1(2),[5 6]) = [PR1' PE1'];
fossilmat2(50:Tmat2(2),[5 6]) = [PR2' PE2'];
fossilmat3(50:Tmat3(2),[5 6]) = [PR3' PE3'];

wholemat0(1:Tmat0(2)+1,[1:3 5 6]) = fossilmat0(:,[1:3 5 6]);
wholemat1(1:Tmat1(2)+1,[1:3 5 6]) = fossilmat1(:,[1:3 5 6]);
wholemat2(1:Tmat2(2)+1,[1:3 5 6]) = fossilmat2(:,[1:3 5 6]);
wholemat3(1:Tmat3(2)+1,[1:3 5 6]) = fossilmat3(:,[1:3 5 6]);

save smthResults.mat
 


